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ABSTRACT 

Adaptive optics (AO) allows one to derive the point spread function (PSF) simultaneously to the science image, 
which is a major advantage in post-processing tasks such as astrometry/photometry or deconvolution. Based on 
the algorithm of Veran et al. (1997), PSF reconstruction has been developed for four different AO systems so far: 
PUEO, ALFA, Lick- AO and Altair. A similar effort is undertaken for NAOS/VLT in a collaboration between 
the group PHASE (Onera and Observatoire de Paris/LESIA) and ESO. In this paper, we first introduce two 
new algorithms that prevent the use of the so-called "Uij functions" to: (1) avoid the storage of a large amount 
| of data (for both new algorithms), (2) shorten the PSF reconstruction computation time (for one of the two) 
■ and (3) provide an estimation of the PSF variability (for the other one). We then identify and explain issues in 
the exploitation of real-time Shack-Hartmann (SH) data for PSF reconstruction, emphasising the large impact 
of thresholding in the accuracy of the phase residual estimation. Finally, we present the data provided by the 
NAOS real-time computer (RTC) to reconstruct PSF {(1) the data presently available, (2) two NAOS software 
modifications that would provide new data to increase the accuracy of the PSF reconstruction and (3) the tests 
Q-f of these modifications) and the PSF reconstruction algorithms we are developing for NAOS on that basis. 
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1. INTRODUCTION 
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The knowledge of the PSF of the instrument is of prime importance to accurately calibrate observations. Nu- 
merous image restoration techniques, such as deconvolution algorithms, make use of the PSF, or its Fourier 
transform, the optical transfert function (OTF), to fully restore the image quality affected by the atmopheric 
turbulence. Astrometric and photometric algorithms, e.g. DAOPHOT (Stetson, 1987) or Starfinder (Diolaiti et 
al., 2000), sometimes coupled to deconvolution algorithms, also need an estimation of the PSF. 

Veran et al. (1997) have demonstrated that AO allows one to derive the PSF from wavefront-related data 
delivered by the AO RTC. They have been the first to develop such a PSF reconstruction algorithm and to 
implement it on an AO system: PUEO, the CFHT curvature AO system (Arsenault et al., 1994). Since then, 
the algorithm is routinely used to provide astronomers with the (on-axis) PSF associated to their observations. 

Based on this algorithm, PSF reconstruction has been developed for three other AO systems, equipped this 
time with SH wavefront sensors (WFS): (1) ALFA (Kasper et al., 2000), the SH AO system of the Calar Alto 
3.5m telescope, by Weiss (2003); (2) Altair (Herriot et al., 2000), the 4-quadrant SH AO system of the Gemini 
North telescope, by Jolissaint et al. (2004); (3) the SH AO system of the UCO/Lick Observatory's 3 m Shane 
Telescope (Bauman et al., 2002), by Fitzgerald (2004). These algorithms have been tested during few runs of 
observations and have lead to good results. 
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With the goal to reconstruct PSF with NAOS, the SH AO system of NACO (Roussct et al., 2000; Lenzen 
et al., 1998) at VLT, a similar effort has been undertaken in a collaboration between the group PHASE (Onera 
and Observatoire de Paris/LESIA) and ESO: PSF reconstruction had indeed been considered since the design 
phase of NAOS. In a first step, we have contacted MPIA (S. Egner and S. Hippler) since they had been the 
first to develop a PSF reconstruction algorithm for a SH AO system. Though, the large differences between 
the two systems (modes of the systems, available wavefront-related measurements) have lead us to elaborate a 
NAOS-dedicated algorithm, the general architecture being always based on the algorithm of Veran et al. (1997). 

In the first section, we introduce new algorithms for PSF reconstruction. In the second section, we identify 
and explain issues in the exploitation of real-time SH data for PSF reconstruction. In the third section, we 
describe the data provided by the NAOS RTC to reconstruct PSF and the contemplated software modifications 
that would provide new data to increase the accuracy of the PSF reconstruction. We also detail the algorithms 
associated with these data. 



2.1. Reminder on the Veran et al. (1997) algorithm 
2.1.1. The long-exposure AO-corrected PSF expression 

Following the PSF reconstruction algorithm developed by Veran et al. (1997), assuming that the phase is quasi- 
stationary over the pupil, the AO-corrected monochromatic long-exposure OTF is decomposed as follows: 



The phase 4> e can be split into two parts: <fr t , which belongs to the vector space spanned by the mirror modes, 
and <f) €± , which is orthogonal to the former space. 



2. NEW ALGORITHMS FOR PSF RECONSTRUCTION 




(1) 



(OTF(p/X)) = (OTF^ (p/X)) x (OTF^ (p/X)) x OTF tel (p/X) 



(2) 



This expression can be written 



(OTF(p/X)) = cxp ( - l -D^ (/J)) x exp ( - (p)) x OTF tel (p/X) 



(3) 



where: 



• (OTF ( f >c (p/X) J is the mean attenuation of the OTF due to the partial correction of AO, 

• OTF tcl (p/X) is the perfect telescope OTF, 

• (OTF^^ (p/A)\ is the mean attenuation of the OTF due the mirror component of the phase, i.e. the 



" mean residual phase OTF" , 
• (OTFfa {p/X) \ is the mean attenuation of the OTF due the component of the phase belonging to the 



space perpendicular to the mirror space, i.e. the "mean perpendicular phase OTF", 
• Z)^ (p) is the mean structure function of the residual phase, i.e. the "mean residual phase structure 



function" , 

• (p) is the mean structure function of the perpendicular phase, 

• p is a pupil plane coordinate vector, 



• A is the wavelength of observation. 



The corresponding AO-corrected monochromatic long-exposure PSF is derived as the Fourier transform of the 
OTF. 



From the expression of the residual phase structure function: 

= ((&„ (2) - (2+ Pi)') (4) 
and the decomposition of the phase on the basis of the mirror modes Mi(x): 

N 

<f> eil (x,t) = J2e lli (t)M i (x) (5) 

one obtains: 

N N 

D^(x,p) = EE<W>( M *( f ) + P))( M M ~Mj(S +p)) (6) 

»=i j=i 

where af a coordinate vector in the pupil plane. 

The mean residual phase structure function (p) is the mean of D<t>, (x, p) over x: 

/ D (x,p)P(x)P(£ + p)dx 

(P) = 1 j (7) 

/ P{x) P{x + p)dx 

where P(x) is the pupil function. 

2.1.2. The Uij(p) functions 
Defining the function Uij (p) as: 

f (Mi(x) - Mi{x + p))(Mj{x) - M 3 (x + p))P(x)P(x + p)dx 

J (8) 



P{x) P{x + p)dx 
Eq. 7, can be rewritten: 

JV N 

%(^) = EEw^w ( 9 ) 
»=i j=i 

This is a key equation for the experimental PSF reconstruction. The covariance matrix (e||e||*) has to be measured 
experimentally on the AO system itself, by averaging the cross-products of wavefront measurements obtained 
during the time of the image exposure. 

In the current PSF reconstruction algorithms, derived from Veran et al. (1997), the matrix (eiieii ) is the 
basic entry point from which one can deduce successively the phase structure function, the OTF, and then the 
PSF. Additionally, one has to compute, store once for all, and also read during the reconstruction process the 
Uij(p) functions. 

In Eqs. 9, the i and j indices play a symmetric role, so that there are actually N x (N + l)/2 "useful" 
Uij(p) functions. As an example, in the case of NAOS, the 159 compensated modes lead to 12720 "useful" Uij(p) 
functions. Today, the large number of Uij(p) hence represents, depending on the array size and data type, several 
gigabytes of data to compute, store and read. This leads to a heavy PSF reconstruction process, which will turn 
out to be impossible to handle in the future since next AO systems are expected to have a largely increased 
number of modes: about 1370 actuators for the VLT-Planet Finder AO system (Fusco et al., 2005), several tens 
of thousands for extremely large telescopes. We propose in the following a way to achieve this computation, 
starting from the same covariance matrix (eiieii*), without using the Uij(p). 



2.2. Theory of the proposed algorithms 

2.2.1. The V it algorithm 

Let's consider the vector e\\* made of the {e\\i}i=i...N coefficients, i.e. the vector representing <j> e ^(x, i) in the 
basis of the mirror modes Mj(x). The eigen decomposition of the residual phase covariance matrix is: 

A = B t (e in t )B (10) 

where A is a diagonal matrix that contains the {Aj}j = i..jv eigenvalues and B is the matrix of eigenvectors: 
B l B = BB l = Id. Equation 10 can be written: 

A=((B t e\\)(B t H ) t ) (11) 

The vector i] equal to B t e\\ represents £|| (x, t) in the basis that diagonalizes the residual phase covariance matrix. 
Its coefficients are noted {?7i}i=i...jv- From Eq. 11, the covariance matrix {ryrf) is diagonal, i.e. in this new basis, 
the residual phase covariance matrix is diagonal and the mean residual phase structure function reduces to: 

JV JV 

i=i i=\ 

where the Vij(p) functions are equivalent in the new basis to the Uij(p) functions (Eq. 8). 

After this change of basis, the computation of the residual phase OTF only requires the computation of a 
number N of functions Vu (p) . Though, these Vu (p) functions have to be computed on the fly for each estimation 
of the mean residual phase structure function. 

2.2.2. The "instantaneous PSF" algorithm 

The solution we propose here is similar to the algorithm presented by Roddier (1990) to simulate atmospherically 
distorted wavefronts, from the covariance matrix of the coefficients of their expansion in Zcrnike modes. We 
extend this algorithm to any modal basis and covariance matrix, and use it to reconstruct AO-corrected PSFs. 

Let's consider again the eigen decomposition of the residual phase covariance matrix: 

(eiieii 4 ) = BAB* (13) 

If one generates a vector i] whose coefficients are independent Gaussian random variables with zero mean and 
variance equal to the eigenvalue A i; i.e. {qif) = A, then the vector (3 = Br] is a set of correlated random variables 
whose covariance matrix is (e||6||*): 

(/3/3*) = {BrirfB 1 ) = BAB* = (e^) (14) 
The phase represented by the vector (3 is: 

JV 

4>{x,t) = Y,mMi{x) (15) 

i=l 

and the "instantaneous" PSF corresponding to that phase is: 

PSF\\(x,t)= ^T(exp(i<73(x,i)) 2 (16) 

Then, by generating random rj vectors such that {rjrf) = A, we build instantaneous PSFs that, in average, 
converge to the long-exposure PSF of the mirror space. Note that the latter is not the "full PSF" as would be 
observed at the telescope since it does not include the uncorrected part of the phase (cf. Eq. 2). 

*to lighten the expressions, e\\(t) will be abbreviated into en 



2.3. Choice between the different algorithms 

The Uij and Vu algorithms mathematically produce exactly the same OTFs. Though, the Vu algorithm requires 
the computation of N functions (where N is the number of modes) whereas the Uij algorithm requires the 
computation of N x (N + l)/2 functions. Even if these N Vu functions have to be computed on the fly, this 
computation is faster than reading N x (N + l)/2 stored Uij functions and it prevents from the storage of a large 
amount of data, so that the Vu algorithm is always to be preferred to the Uij one. 

As noticed by Conan (1994), averaging short-exposure OTFs, as we do in the "instantaneous PSF" algorithm, 
is a process that converges very slowly, especially at large D/ro or low correction level. In addition, it does not 
lead to the infinitely long exposure OTF since a given number of short exposures are averaged. Besides, Conan 
(1994) has shown that in the poor correction case, the error in computing the long-exposure OTF of such an 
algorithm is larger than for the Uij algorithm, and consequently the Vu algorithm as well. 

Though, we emphasize that, in addition to the OTF computation itself, the "instantaneous PSF" algorithm 
can provide an estimation about the variability of the OTF, that can help a lot in some deconvolution algorithms. 
The estimation of the infinitely long exposure OTF that results from the convergence of our "instantaneous PSF" 
algorithm and that corresponds to a given covariance matrix (e||£||*) is unique: let us call it OTF^p/X). The 
dispersion in the random generation of OTFs can be computed as 



where OTFi is the i th draw of a randomly-generated OTF. If we call OTF ohs (p/\) the OTF actually observed 
on the instrument during the given, non infinite, observing time T; nt , when the given covariance matrix (e||e||*) 
was measured, we can evaluate how far our estimation OTF^ is from OTF Q ^ s by writing: 



where n is the equivalent number of independent realisations of PSFs, whose sum has resulted in the final PSF 
observed by the instrument during the given, non infinite, integration time Tj nt . An estimation of n could be 
obtained for example from a full simulation of the AO system under the same atmospheric conditions as during 
the observation. It can also be reasonnable to consider that the impact of the correction by the AO system is 
to shorten the image correlation time compared to the image correlation time To (A) of the atmospheric seeing 
(Rigaut et al., 1991). We can then find a lower bound given by n > Ti nt /r (A). 



3. ISSUES IN THE EXPLOITATION OF REAL-TIME SHACK-HARTMANN DATA 



3.1. Description of the problem 
3.1.1. Using the Yao simulation software 

In order to test the algorithms introduced in Sect. 2.2, we have used the end-to-end AO simulation software Yao 
(version 3.6), written by Rigaut (2005). We have tested our algorithms in the simple case of a 6x6 subpupil 
SH by using the default configuration file " sh6m2-bench.par" delivered in the Yao installation package. The 
correction was performed with a 45 actuator piezostack deformable mirror (DM) plus a tip-tilt mirror. The 
correction and observation wavelengths of the simulation were 0.65 and 1.65 /im respectively. D/ro was equal to 
42.44 at 0.5 fim. 

Since we aimed at testing our algorithms with different conditions of correction, we have run the simulation 
with a guide star magnitude ranging from 5 to 14 so that the resulting Strehl ratio was ranging from 0.54 
down to 0.01 (Fig. 1 left). We have stored all the values e\\(t) obtained from the simulation, computed the 
covariance matrix (e||e||*) and ran the Uij, Vu and "instantaneous PSF" algorithms to derive the corresponding 
"atmospheric" OTFs (i.e. not multiplied by the telescope OTF; Fig. 1 right). 

The Strehl ratios that one would expect from the OTF profiles (Fig. 1 right), i.e. computed from the WFS 
data, do not correspond to the Strehl ratios computed from the PSFs resulting from the simulation (Fig. 1 left), 




(17) 




(18) 



FOR PSF RECONSTRUCTION 




Figure 1. Left: Strehl ratio at 1.65 /im vs. guide star V magnitude, from Yao. Yao computes the Strehl ratios from 
the PSFs resulting from the simulation (as would be observed on a camera). Right: Circular mean of the "atmospheric" 
OTFs for guide star magnitudes of 5, 12 and 14 (from top to bottom), computed from the WFS data. At this scale, the 
profiles obtained with the different algorithms (Uij, Vu and "instantaneous PSF") are almost superimposed. 

even if one takes into account that these OTFs only correspond to the residual phase, without correction for the 
noise, the aliasing and the perpendicular parts of the phase. This is particularly evident for the 14 th guide star 
magnitude case, for which one would expect much smaller OTF values for a large range of normalized frequencies. 
In order to understand this discrepancy, we have reproduced the test with another AO simulation software. 

3.1.2. Using the ONERA AO simulation software 

To trace the origin of the problem, we have tested the algorithm with an end-to-end Monte Carlo-based AO 
simulation software developed at ONERA (Conan et al., 2004). We have tested our algorithms in the simple 
case of a NAOS-like AO system (Rousset et al., 2000): a 14x14 subpupil SH with 8x8 pixels per sub-aperture, 
a read-out noise of 3 e~ per pixel and a sampling frequency of 500 Hz. The correction was performed with a 
185 actuator piezostack DM plus a tip-tilt mirror. The correction and observation wavelengths of the simulation 
were 0.65 and 2.2 ^m respectively. The seeing was 0.85" at 0.5 /im. 

Similarly to what we have done with Yao, we have run the simulation with a guide star magnitude ranging 
from 7 to 15, leading to a resulting Strehl ratio from ~70% down to ^0.2% (Fig. 2 left), and for a given guide 
star magnitude, we have stored the values eii(t) obtained from the simulation, computed the covariance matrix 
(e||e||*) and run the Uij, Vu and "instantaneous PSF" algorithms to derive the corresponding "atmospheric" 
OTFs (Fig. 2 right). 

The behaviour we had encountered with Yao has still been observed in this other simulation: there exists a 
discrepancy between the Strehl ratios computed from the "imaging path" -related data (the PSFs, left graph of 
Fig. 2) and the Strehl ratios expected from the WFS data (the OTFs, right graph of Fig. 2). 

3.2. Origin of the problem 

We have identified unambiguously that the unexpected high OTF values are related to the correlation between 
the error and the turbulence signal. 

In order to analyse the effect in more detail, we have simulated the behaviour of the SH measurement in 
one subaperture in particular. We simulate the incoming wavefront, compute the associated image with a high 
resolution and derive the position of its centre of gravity (hereafter the "true centroid position"). Then, we 
simulate the technical characteristics of the sensor, namely: the averaging of the image over the pixel area, the 
limited number of pixels and field of view, the detector noise, the signal photon noise. We apply a threshold on 
the image and compute the centre of gravity (hereafter the "measured centroid position"). After comparing the 
measured centroid position (s m , for measured slope) with the true one (s ), it appears clearly that: 

• the error s m — s is highly correlated to the signal value (i.e. to the true centroid position), which is 
contradictory with the assumptions usually made in the literature; 

• the correlation is, among others, a function of the threshold level. 
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Figure 2. Left: Strehl ratio at 2.2 /im vs. guide star visual magnitude for a NAOS-like AO system. As for Yao, 
the Strehl ratios are computed from the PSFs obtained from the simulation (as would be seen in the "imaging path"). 
Right: Circular mean of the "atmospheric" residual OTFs for guide star magnitudes of 7.3, 12.3 and 13.6 (from top to 
bottom), computed from the WFS data. At this scale, the profiles obtained with the different algorithms (Uij, Vu and 
"instantaneous PSF") are almost superimposed. 



Therefore, instead of the usual equation s„ 
written under the form: 



so + n, the measured signal in a subaperture should rather be 
s m = G.sq + n (19) 



where G is a coefficient and n an additive noise, uncorrelated with the signal sq- The value of G can be found 
with a linear regression (least-squares fit by a function ax + b) of the measured data s m with respect to so- 
The fit residuals (i.e. n) are, de facto, uncorrelated with so, as a consequence of the least-squares approach. In 
addition, it should be noticed that n may not have a zero mean, depending on the location of the reference slope 
in the subaperture. The value of G is related to the noisy pixel values that peak above the threshold and bias 
the estimation of the centroid location, because they weight in the denominator of the expression of the centre 
of gravity, given by (in the x- measurement case): 
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(20) 



where Si, iVj and Xi are the flux signal, the flux noise and the position of the pixel i in the given supaperture. 
Equation 20 is comparable to Eq. 19, where the mean value of G can be expressed in terms of the total flux 
Nphot detected above the threshold, and in terms of the sum of the noise N noise above the threshold level 

G~ N N Tn < 21 > 

lyj phot T" J* noise 



We have plotted in Fig. 3 the value of G versus the threshold level. While Eq. 21 can closely predict the 
value of G for the low threshold values, it cannot predict the behavior for higher threshold levels, where the 
simulations show that, in some cases, we may even have G > 1. The interesting point is that the error reaches 
a minimum when the correlation is zero, which coincides with a regression coefficient of G = 1. 

We have reproduced this particular " gain effect" , encountered in our end-to-end simulations described in 
Sect. 3.1, using the ONERA AO simulation software, not in its end-to-end configuration but: 

• the slope is computed as a simple phase difference at the edge of the corresponding subpupil, 

• we then multiply by a gain factor G = Nphot I (Nphot + N noise ), where N noise follows a truncated Gaussian 
statistics in order to reproduce the threshold effect. 

• the noise is then added to the slope measurement. 
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Figure 3. Left : Coefficient G vs. threshold from the linear regression of the measured data s m through G.s + constant. 
Middle : Correlation coefficient between the signal so and the error s m — so vs. threshold. Right : Variance of the error 
((s m — so) 2 ) vs. threshold. Conditions for the 3 figures: 7x7 subapertures, 8-m telescope, r<j = f3 cm at 0.6 microns, 
detector read-out noise: 3 e~, pixel scale: 0.4" /pixel, centroid computed on 6 x 6 pixels, sampling frequency: 50 Hz, star 
magnitude: R=I 5, AA = 300 nm, global throughput: 0.25. 

Indeed, the resulting upper dashed curve on the left graph of Fig. 4 reproduces the behaviour of our two end- 
to-end full-propagation simulations: the Strehl ratio reaches a lower limit at low flux instead of getting down 
to zero as observed from the "imaging path" -related data (lower continuous curve on the left graph of Fig. 4). 
However, if we correct the slope measurements from the gain of the centroid measurement, we find the expected 
behaviour (central dotted curve on the left graph of Fig. 4, which is not superimposed on the lower continuous 
curve because the slope measurements have not been corrected for the noise). 

As a conclusion, the threshold level is a key issue when reconstructing PSFs from real SH data. An improper 
threshold level may bias the data, leading to a wrong estimation of the noise and of the phase residuals. This issue 
may also impact on all the tasks that make use of the slope measurements, such as performance estimations (as 
done for NAOS), modal optimisation, etc. Surprisingly, its impact on the performance of the AO compensation 
itself is maybe less sensitive, since the gain G acts exactly as the close-loop gain would do: it directly controls 
the system bandwidth, and only makes the system "slower" when G < 1. 

4. NAOS DATA FOR PSF RECONSTRUCTION AND ASSOCIATED ALGORITHMS 
4.1. Data presently available 

The delivery of wavefront-related data to estimate the PSF has been among the NAOS main specifications to 
maximise the scientific returns of the instrument. In that purpose, the NAOS RTC computes over the image 
integration time the mean of the following data (Rabaud et al., 2000): 

• at WFS frequency: Zernike Mode Mean, Zernike Mode Autocorrelation (128 first points, mean value not 
subtracted), Intensity Mean, Intensity Variance, Image Mean, Image Variance, Image pixel null/no-null 
count; 

• at 50 Hz: mean of the modal coefficients deduced from the residual slopes (e, hereafter residual modal 
mean), covariance of the modal coefficients deduced from the residual slopes (C ££ , hereafter residual modal 
covariance matrix) , mean of the modal coefficients deduced from the voltages and covariance of the modal 
coefficients deduced from the voltages (hereafter mirror modal mean and covariance matrix resp.). 

In practice, together with the CONICA FITS images, the observer gets: 

• attached to the images: the two residual modal covariance and mirror modal covariance matrices and the 
corresponding two means; 

• written in the FITS header images: turbulence parameters, such as r , L (cf. Fusco et al. 2004 for details 
on their estimation), and AO loop parameters such as the Zernike mean noise n\. 
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Figure 4. Left: Strehl ratio vs. guide star flux (in number of photons). The slope is here computed as a simple phase 
difference at the subpupil edges and (1) we reproduce the threshold effect (upper dashed curve) and (2) we correct the 
slope measurements from the gain of the centroid measurement, but not from the noise (central dotted curve). These 
Strehl ratios are computed from WFS data for these curves and from "imaging path" -related data for the lower continuous 
one. Right: Same as Fig. 2 right, after correcting the slope measurements from the gain of the centroid measurement. 



4.2. From the available data to the PSF 

4.2.1. Computation of the mean residual phase OTF 

The WFS measurement, ey can be decomposed into: 



(22) 



where en, n and r are the vectors representing the residual phase, the noise and the aliasing in the basis of the 
miror modes. Following Veran et al. (1997), under the assumption of high temporal bandwidth, Eq. 22 leads to: 



(e\ l e ll t ) = (e ll e\)-(nn t ) + (rr t ) 



(23) 



- is directly obtained from the NAOS data: = C ee — (ee*) 

- (rr*) is in a first step neglected. In the future, it will be computed by using a SH simulation of NAOS 

- (nn*) was contemplated to be computed from the Zernike mean noise value nj, given in the CONICA image 
FITS header. To obtain this value, the vector of the residual phase is expressed in the Zernike mode basis, 
the open-loop measurement vector is reconstructed, the noise variances for each Zernike are calculated with the 
auto-correlation method (Gendron & Lena, 1995) and then summed and square-rooted. To take advantage of 
this Zernike mean noise value, one has to assume that the noise is not correlated from a subpupil to another and 
that the Zernike mean noise is equally distributed over all the considered Zernike modes (!). In that case, (nn*) 
is given by (hereafter "n^ algorithm"): 

(nn') = — \ M SM M% M (24) 



2 

SZi. 



where M$m is the "slope to mirror mode" transformation matrix, Mgz the "slope to Zernike mode" transfor- 
mation matrix and M$zi,j its element along line i and column j. 

4.2.2. Computation of the mean perpendicular phase OTF 

As proposed by Veran et al. (1997), the mean perpendicular phase OTF can be computed from simulated phases 
screens, computed at D/r — 1, of which one has subtracted their mirror components to obtain <\> tl _ (x, t). These 
"perpendicular phase screens" have then to be scaled at the proper D/ro value by multiplying them by (D/ro) 5 ^ 6 - 
The ro value of the observations can be directly read in the CONICA image FITS header. Alternatively, it could 
be derived from the mirror modal mean and covariance matrix attached to the CONICA image FITS file. The 
former solution is presently adopted. 

From these scaled perpendicular phase screens, one can derive the mean perpendicular phase OTF with two 
different algorithms: 



a "PSF-like" algorithm: one computes for each perpendicular phase screen the corresponding "instan- 

2 



taneous perpendicular PSF" by 



TTT{^ exp(z0 £i (x, t)) , averages them to get PSF±(x,t) and, after a 



Fourier transform, the perpendicular OTF multiplied by the telescope OTF; 

an " t/ij-like" algorithm: similarly to the commonly used computation of the Uij functions, 



P(x)P(x + p)D^ ± {x,p)dx 

D^ { {x,p) = J - (25) 

/ P(x)P(x + p)dx 

f P(x)P(x + p) (<£ 2 ± (x , t) + <t> 2 e± (x + p,t)- 2<j> e± (x, t)<f> e± (x + p, t)) dx\ 

(26) 



P(x)P(x + p)dx 

(Cor{cj > 2 ± P,P)+Cor{P^ 2 ^P) - 2Cor(^P,^ x P)) 
Cor(P, P) 

T- 1 ({T{^ ± P)T*{P)+HP)^(4>l ± P) - 2F(<j> e± P)F(<i> e± P) 



t- 1 (| HP) I 2 ) 

T- 1 ( y 2(K(H4>l ± P)F*(P))-\T{cj > ^P) | 2 ) 
F- 1 (I W) I 2 ) 



(27) 
(28) 

(29) 



where Cor(f,g) = J f(x).g(x + p)dx, T is the Fourier transform, T 1 the inverse Fourier transform, 5R 
the real part of a complex value and * its conjugate. 

4.3. Contemplated NAOS software modifications 
4.3.1. Description 

The Zernike mean noise value n 2 z is the only NAOS data available to compute the noise part of the OTF 
(Sect. 4.1). Though, as explained in Sect. 4.2.1, the assumptions made in this computation are fairly important 
and might lead to large uncertainties in the PSF estimation. It has lead us to propose to ESO two different RTC 
software modifications with different levels of impact on the software: 

1. providing the vector of variance noises n 2 z . for all considered Zernikes, instead of their "mean" value n 2 
(hereafter "n 2 .. algorithm"). This is the modification that has the slightest impact on the software: since 
the vector of variance noise for all considered Zernikes is already computed in the software to derive n 2 ., 
this software modification would only require to output the considered vector, in the CONICA image FITS 
header for example, instead of outputting n 2 . Always assuming that the noise is not correlated from a 
subpupil to another, we would first solve the equations n 2 . = ^ Mg Z . ,n 2 Zj where n 2 . is the i^-element of 

3 

the vector of variance noises for each slope mesurements. This can be done with a dedicated procedure, e.g. 
the IDL "constrained_min.pro" procedure, (nn*) is then derived using the M$m "slope to mirror mode" 
transformation matrix, 

2. providing the vector of variance noises nf for all considered mirror modes, i.e. the diagonal of (nn*) 
(hereafter "n 2 algorithm"). In this larger software modification scheme, always assuming no correlation 
between the subpupils, we would first solve the equations n 2 . = ^ Mg M . .n 2 , similarly to what is proposed 

j 

for the first software modification, and derive the entire (nn 4 ) covariance matrix with M$m- 
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Figure 5. Drawn (black) and estimated (red) n 2 . slope variance noises with the nl (left), n 2 . (middle) and n 2 (right) 
algorithms, for a given statistics distribution. 



4.3.2. Test of the modifications 

To test the contemplated modifications, we have drawn vectors of 144 n 2 . slope variance noises with several 
statistics distribution, taking into account that the edge subpupils are more noisy than the central ones. Assuming 
no correlation between subpupils, i.e. (n s n* ) is diagonal and made of the n s . elements, we have then: 

• computed (nn ) with the "slope to mode" transformation matrix, and then derived with the Vu algorithm 
the corresponding OTF (hereafter "reference OTF"), that we estimate with the n 2 , n 2 . and n 2 algorithms; 

• computed n 2 with the "slope to Zernike" transformation matrix, used the "n 2 algorithm" to get the 
estimated n 2 . slope variance noises (Fig. 5 left), estimated (nn 1 ) with the "slope to mode" transformation 
matrix, and then derived with the Va algorithm the corresponding estimated OTF; 

• computed nf. with the "slope to Zernike" transformation matrix, used the "n 2 . algorithm" to get the esti- 
mated n 2 . slope variance noises (Fig. 5 middle), estimated (nn*) with the "slope to mode" transformation 
matrix, and then derived with the Vu algorithm the corresponding estimated OTF; 

• computed n 2 with the " slope to mode" transformation matrix, used the " nf algorithm" to get the estimated 
n 2 . slope variance noises (Fig. 5 right), estimated (nn*) with the "slope to mode" transformation matrix, 
and then derived with the Vu algorithm the corresponding estimated OTF; 

To measure the potential gain of the contemplated modifications, we have represented in Fig. 6 the histogram, 
in number of pixels in the OTF, of the relative error between the "reference OTF" and the "estimated OTFs", 
for the statistics distribution of the slope variance noises shown in Fig. 5. We thus show the substantial gain of 
the simplest software modification (providing the vector of variance noises n 2 .) compared to what is presently 
delivered by the NAOS RTC (n 2 ). The gain of the second contemplated modification (providing the vector of 
variance noises n 2 ) seems to be poor with respect to the first one. 
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Figure 6. Histogram of the relative error between the "reference OTF" and the estimated OTF estimated with the rcj 
(upper red curve), n 2 . (middle black curve) and n 2 (lower blue curve) algorithms, for the statistics distribution of the 
slope variance noises given in Fig. 5. 



5. FUTURE WORK 

Regarding the new introduced algorithms, accurate AO simulations arc needed to provide an estimation of n, the 
equivalent number of independent realisations of PSFs whose sum has resulted in the observed PSF. The derived 
OTF variability would then be used in deconvolution algorithm, similarly to what has been done with PUEO 
and MISTRAL (Fusco et al., 1999). Concerning the NAOS PSF reconstruction software itself, the decision to 
proceed with the slightest software modification is still in discussion with ESO. Though, tests of the software with 
the present NAOS RTC data have just begun and more will be accumulated before a release to the community. 
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